#include <iostream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>
#include "shape_align.h"

using namespace std;


class AlignmentKernelStatus
{
  float fMscore;
  float fRscore;
  float fCscore;
  bool fMb;
  bool fRb;
  bool fCb;
  float rMscore;
  float rRscore;
  float rCscore;
  bool rMb;
  bool rRb;
  bool rCb;

  int RowNumGap;
  int ColumnNumGap;

public:
  AlignmentKernelStatus():
    fMb(false),fRb(false),fCb(false),rMb(false),rRb(false),rCb(false){}
  ~AlignmentKernelStatus(){}

  void setfMscore(float ifMscore){fMscore=ifMscore; fMb=true;}
  void setfRscore(float ifRscore){fRscore=ifRscore; fRb=true;}
  void setfCscore(float ifCscore){fCscore=ifCscore; fCb=true;}
  void setrMscore(float irMscore){rMscore=irMscore; rMb=true;}
  void setrRscore(float irRscore){rRscore=irRscore; rRb=true;}
  void setrCscore(float irCscore){rCscore=irCscore; rCb=true;}
  void setRowNumGap(int iRowNumGap){RowNumGap=iRowNumGap;}
  void setColumnNumGap(int iColumnNumGap){ColumnNumGap=iColumnNumGap;}

  float getfMscore() const {return fMscore;}
  float getfRscore() const {return fRscore;}
  float getfCscore() const {return fCscore;}
  float getrMscore() const {return rMscore;}
  float getrRscore() const {return rRscore;}
  float getrCscore() const {return rCscore;}
  bool getfMb() const {return fMb;}
  bool getfRb() const {return fRb;}
  bool getfCb() const {return fCb;}
  bool getrMb() const {return rMb;}
  bool getrRb() const {return rRb;}
  bool getrCb()const {return rCb;}

  vector<float> getfscores() const
  {
    vector<float> fscores;
    if(fMb == true){ fscores.push_back(fMscore); }
    if(fRb == true){ fscores.push_back(fRscore); }
    if(fCb == true){ fscores.push_back(fCscore); }
    if(fscores.size() == 0){ throw "getfscores error"; /*cerr << "getfscores error" << endl; exit(1);*/}
    sort(fscores.begin(),fscores.end(),std::greater<float>());
    return fscores;
  }
  vector<float> getrscores() const
  {
    vector<float> rscores;
    if(rMb == true){ rscores.push_back(fMscore); }
    if(rRb == true){ rscores.push_back(fRscore); }
    if(rCb == true){ rscores.push_back(fCscore); }
    if(rscores.size() == 0){throw "getfscores error"; /*cerr << "getfscores error" << endl; exit(1);*/}
    sort(rscores.begin(),rscores.end(),std::greater<float>());
    return rscores;
  }

  int getRowNumGap() const {return RowNumGap;}
  int getColumnNumGap() const {return ColumnNumGap;}
  int getNumGap() const {return RowNumGap+ColumnNumGap;}
};

class PathStatus
{
  double RowCoverage;
  double ColumnCoverage;
  double SumDelta;
  int BackDirection; //Diagonal:0, RowMove:1, ColumnMove:2 
  int RowNumGap;
  int ColumnNumGap;

public:

  ~PathStatus(){}
  PathStatus() : SumDelta(100000000),RowNumGap(0),ColumnNumGap(0) {}
  void setRowCoverage(double iRowCoverage){RowCoverage=iRowCoverage;}
  void setColumnCoverage(double iColumnCoverage){ColumnCoverage=iColumnCoverage;}
  void setSumDelta(double iSumDelta){SumDelta=iSumDelta;}
  void setBackDirection(int iBackDirection){BackDirection=iBackDirection;}
  void setRowNumGap(int iRowNumGap){RowNumGap=iRowNumGap;}
  void setColumnNumGap(int iColumnNumGap){ColumnNumGap=iColumnNumGap;}
  double getRowCoverage() const {return RowCoverage;}
  double getColumnCoverage() const {return ColumnCoverage;}
  double getSumDelta() const{return SumDelta;}
  int getBackDirection() const{return BackDirection;}
  int getRowNumGap() const {return RowNumGap;}
  int getColumnNumGap() const {return ColumnNumGap;}
  int getNumGap() const {return RowNumGap+ColumnNumGap; }
};
/*sortの基準設定*/
bool operator<( const PathStatus &PS1, const PathStatus &PS2) { return PS1.getSumDelta() > PS2.getSumDelta(); } //これで降順になる??        
bool operator==( const PathStatus &PS1, const PathStatus &PS2) { return PS1.getSumDelta() == PS2.getSumDelta(); }
/* --- */

template <class Status> //Status = PathStatus or AlignmentKernelStatus
class CELL{

  vector<Status> SS;
  vector<bool> bools;
  vector<int> Index;

public:

  ~CELL(){}
  CELL(int heigth) : SS(heigth),bools(heigth+1,false),Index(0) {}
  void setIndex(int i){Index.push_back(i);}
  void setTrue(int i){bools.at(i)=true;}
  vector<Status> &getSS() {return SS;}
  const vector<Status> &getSS() const {return SS;}
  const vector<bool> &getbools() const {return bools;}
  bool getbool(int i) const {return bools.at(i);}
  const vector<int> &getIndex() const {return Index;}

  void setTrueAndIndex(int i)
  {
    if(this->getbools().at(i) == false)
      {
	this->setTrue(i);
	this->setIndex(i);
      }
  }
};

class DPs{

protected:

  const vector<double> &RowCoverages;
  const vector<double> &ColumnCoverages;
  int MaxRowGap;
  int MaxColumnGap;
  bool BothEndZeroOption;

  double makeRowGapCoverage(int Row) const;
  double makeColumnGapCoverage(int Col) const;
  double GapScoreOfRowMove(int Row,int Col) const;
  double GapScoreOfColumnMove(int Row, int Col) const;
  double GapScoreOfDiagonalMove(int Row, int Col) const;

  // use when Kernel Alignment
  double beta;
  float lnSum(const vector<float>& scores) const;
  float makefMscore(const vector<float>& scores, int Row, int Col) const;
  float makefRscore(const vector<float>& scores, int Row, int Col) const;
  float makefCscore(const vector<float>& scores, int Row, int Col) const;
  float makerscore(const vector<float> &scores) const;
  // 

public:
  ~DPs(){;}
  DPs(const vector<double>& iRowCoverages, const vector<double>& iColumnCoverages, int iMaxRowGap, double ibeta=-1.0, bool iBothEndZeroOption=false) 
    : RowCoverages(iRowCoverages),ColumnCoverages(iColumnCoverages), MaxRowGap(iMaxRowGap), beta(ibeta),
      BothEndZeroOption(iBothEndZeroOption)
  {
    int RowSize = iRowCoverages.size()+1;
    int ColumnSize = iColumnCoverages.size()+1;
    int MaxGap = abs( RowSize-ColumnSize );
    // if(GapOption == "flow")
    //   {
    //     MaxRowGap = (int) iColumnCoverages.size()*alpha;
    //     if(MaxRowGap >= RowSize){MaxRowGap=RowSize-1;}
    //   }
    // else if(GapOption == "fix")
    //   {
    //     MaxRowGap = fixGap; 
    //     if(MaxRowGap >= RowSize){MaxRowGap=RowSize-1;}
    //   }
    if(MaxRowGap >= RowSize){MaxRowGap=RowSize-1;}

    MaxColumnGap = MaxGap+MaxRowGap;
    if(MaxColumnGap >= ColumnSize){MaxColumnGap=ColumnSize-1;}
  }
};

class ShapeAlignment1 : public DPs
{
  vector< vector< CELL<PathStatus> > > cells;

  void DiagonalMove(int Row,int Col);
  void RowMove(int Row,int Col);
  void ColumnMove(int Row,int Col);

public:

  ~ShapeAlignment1(){}
  ShapeAlignment1(const vector<double>& iRowCoverages, const vector<double>& iColumnCoverages, int iMaxRowGap)
    : DPs(iRowCoverages,iColumnCoverages, iMaxRowGap)
  {
    int RowSize = iRowCoverages.size()+1;
    int ColumnSize = iColumnCoverages.size()+1;
    cells.assign( RowSize, vector< CELL<PathStatus> >(ColumnSize, CELL<PathStatus>(MaxRowGap+MaxColumnGap+1) ) ); 
  }
  
  void init();
  void Alignment();
  void TraceBack(vector<double> &Alignment_RowVector, vector<double> &Alignment_ColumnVector);

};

class KernelAlignment1 : public DPs
{
  vector<vector< CELL<AlignmentKernelStatus> > > fcells; //use in forward
  vector<vector< CELL<AlignmentKernelStatus> > > rcells; //use in Backward
  
  void initForwardRow();
  void initForwardCol();
  void CalfMscores(int Row, int Col);
  void CalfRscores(int Row, int Col);
  void CalfCscores(int Row, int Col);

  void initBackwardRow();
  void initBackwardCol();
  void Calrscores(int Row, int Col);
  void setprerMscoreInVector(vector<float> &prerscores, int Row, int Col, int index);
  void setprerRscoreInVector(vector<float> &prerscores, int Row, int Col, int index);
  void setprerCscoreInVector(vector<float> &prerscores, int Row, int Col, int index);

  float AllProb;
  vector< vector<float> > MpostProbs;
  void CalAllProb();
  float CalMpostProb(int Row, int Col);

public:

  ~KernelAlignment1(){}
  KernelAlignment1(const vector<double>& iRowCoverages, const vector<double>& iColumnCoverages, int iMaxRowGap, double ibeta=-1.0, bool iBothEndZeroOption=false)
    : DPs(iRowCoverages,iColumnCoverages, iMaxRowGap, ibeta, iBothEndZeroOption)
  {
    int RowSize = iRowCoverages.size()+1;
    int ColumnSize = iColumnCoverages.size()+1;
    fcells.assign( RowSize, vector< CELL<AlignmentKernelStatus> >(ColumnSize, CELL<AlignmentKernelStatus>(MaxRowGap+MaxColumnGap+1) ) ); 
    rcells.assign( RowSize, vector< CELL<AlignmentKernelStatus> >(ColumnSize, CELL<AlignmentKernelStatus>(MaxRowGap+MaxColumnGap+1) ) ); 
    MpostProbs.assign(iRowCoverages.size(), vector<float>(iColumnCoverages.size(), 0) );
  }

  void initForward();
  void Forward();
  void initBackward();
  void Backward();
  void CalMPostProbs();
  vector< vector<float> > getMpostProbs() const {return MpostProbs;}  
};

class KernelAlignment2 : public DPs // use when alpha=1.0
{
  vector<vector< CELL<AlignmentKernelStatus> > > fcells; //use in forward
  vector<vector< CELL<AlignmentKernelStatus> > > rcells; //use in Backward
  
  void initForwardRow();
  void initForwardCol();
  void CalfMscores(int Row, int Col);
  void CalfRscores(int Row, int Col);
  void CalfCscores(int Row, int Col);

  void initBackwardRow();
  void initBackwardCol();
  void Calrscores(int Row, int Col);
  void setprerMscoreInVector(vector<float> &prerscores, int Row, int Col);
  void setprerRscoreInVector(vector<float> &prerscores, int Row, int Col);
  void setprerCscoreInVector(vector<float> &prerscores, int Row, int Col);

  float AllProb;
  vector< vector<float> > MpostProbs;
  void CalAllProb();
  float CalMpostProb(int Row, int Col);

public:

  ~KernelAlignment2(){}
  KernelAlignment2(const vector<double>& iRowCoverages, const vector<double>& iColumnCoverages, int iMaxRowGap, double ibeta=-1.0, bool iBothEndZeroOption=false)
    : DPs(iRowCoverages,iColumnCoverages, iMaxRowGap, ibeta, iBothEndZeroOption)
  {
    int RowSize = iRowCoverages.size()+1;
    int ColumnSize = iColumnCoverages.size()+1;
    fcells.assign( RowSize, vector< CELL<AlignmentKernelStatus> >(ColumnSize, CELL<AlignmentKernelStatus>(1) ) ); 
    rcells.assign( RowSize, vector< CELL<AlignmentKernelStatus> >(ColumnSize, CELL<AlignmentKernelStatus>(1) ) ); 
    MpostProbs.assign(iRowCoverages.size(), vector<float>(iColumnCoverages.size(), 0) );
  }

  void initForward();
  void Forward();
  void initBackward();
  void Backward();
  void CalMPostProbs();
  vector< vector<float> > getMpostProbs() const {return MpostProbs;}  
};

double DPs::makeRowGapCoverage(int Row) const
{
  if(Row==0)
    {
      if(BothEndZeroOption){return 0;}
      return RowCoverages.at(0);
    }
  else if(Row==RowCoverages.size())
    {
      if(BothEndZeroOption){return 0;}
      return RowCoverages.back();
    }
  else{return ( RowCoverages.at(Row) + RowCoverages.at(Row-1) )/2;}
}

double DPs::makeColumnGapCoverage(int Col) const
{
  if(Col==0)
    {
      if(BothEndZeroOption){return 0;}
      return ColumnCoverages.at(0);
    }
  else if(Col==ColumnCoverages.size())
    {
      if(BothEndZeroOption){return 0;}
      return ColumnCoverages.back();
    }
  else{return ( ColumnCoverages.at(Col) + ColumnCoverages.at(Col-1) )/2;}
}

double DPs::GapScoreOfRowMove(int Row,int Col) const
{
  return abs( RowCoverages.at(Row-1)-makeColumnGapCoverage(Col) ); 
}

double DPs::GapScoreOfColumnMove(int Row, int Col) const
{
  return abs( makeRowGapCoverage(Row)-ColumnCoverages.at(Col-1) ); 
}

double DPs::GapScoreOfDiagonalMove(int Row, int Col) const
{
  return abs( RowCoverages.at(Row-1)-ColumnCoverages.at(Col-1) );
}

void ShapeAlignment1::DiagonalMove(int Row,int Col)
{
  CELL<PathStatus> &precell=cells.at(Row-1).at(Col-1);
  CELL<PathStatus> &Nowcell=cells.at(Row).at(Col);
  
  vector<int>::const_iterator vit = precell.getIndex().begin();
  while( vit != precell.getIndex().end() ) 
    {
      const PathStatus &prePS=precell.getSS().at(*vit);
      int NumGap=prePS.getRowNumGap()+prePS.getColumnNumGap();
      double iSumDelta=prePS.getSumDelta()+GapScoreOfDiagonalMove(Row,Col);
      PathStatus &NowPS=Nowcell.getSS().at(NumGap);		 

      if(iSumDelta <= NowPS.getSumDelta())
	{
	  NowPS.setSumDelta( iSumDelta );
	  NowPS.setRowCoverage( RowCoverages.at(Row-1) );
	  NowPS.setColumnCoverage( ColumnCoverages.at(Col-1) );
	  NowPS.setBackDirection( 0 );
	  NowPS.setRowNumGap( prePS.getRowNumGap() );
	  NowPS.setColumnNumGap( prePS.getColumnNumGap() );
	  
	  Nowcell.setTrueAndIndex(NumGap);
	}
      ++vit;
    }
}

void ShapeAlignment1::RowMove(int Row,int Col)
{
  CELL<PathStatus> &precell=cells.at(Row-1).at(Col);
  //if(Col==0){cerr << Row << " : " << Col << " : " << cells.size() << " : " << MaxRowGap << " : " << MaxColumnGap << endl;}  
  CELL<PathStatus> &Nowcell=cells.at(Row).at(Col);

  vector<int>::const_iterator vit = precell.getIndex().begin();
  while( vit != precell.getIndex().end() ) 
    {
      PathStatus &prePS=precell.getSS().at(*vit);
      int precellRowGap=prePS.getRowNumGap();

      if( precellRowGap < MaxRowGap ) 
	{
	  double iSumDelta=prePS.getSumDelta()+GapScoreOfRowMove(Row,Col);
	  int NumGap=prePS.getRowNumGap()+prePS.getColumnNumGap()+1;

	  //if(Col==0)cerr << MaxRowGap << " : " << MaxColumnGap << " : " << NumGap << " : " << iSumDelta << endl;
	  PathStatus &NowPS=Nowcell.getSS().at(NumGap);		

	  if(iSumDelta <= NowPS.getSumDelta())
	    {
	      NowPS.setSumDelta( iSumDelta );
	      NowPS.setRowCoverage( RowCoverages.at(Row-1) );
	      NowPS.setColumnCoverage( makeColumnGapCoverage(Col) );
	      NowPS.setBackDirection( 1 );
	      NowPS.setRowNumGap( prePS.getRowNumGap()+1 );
	      NowPS.setColumnNumGap( prePS.getColumnNumGap() );
	      
	      Nowcell.setTrueAndIndex(NumGap);
	    }
	}
      ++vit;
    }
}

void ShapeAlignment1::ColumnMove(int Row,int Col)
{
  CELL<PathStatus> &precell=cells.at(Row).at(Col-1);
  CELL<PathStatus> &Nowcell=cells.at(Row).at(Col);
  
  vector<int>::const_iterator vit = precell.getIndex().begin();
  while( vit != precell.getIndex().end() ) 
    {
      PathStatus &prePS=precell.getSS().at(*vit);
      int precellColumnGap=prePS.getColumnNumGap();
      if( precellColumnGap < MaxColumnGap ) 
	{
	  double iSumDelta=prePS.getSumDelta()+GapScoreOfColumnMove(Row,Col);
	  int NumGap=prePS.getRowNumGap()+prePS.getColumnNumGap()+1;
	  PathStatus &NowPS=Nowcell.getSS().at(NumGap);		 
	  
	  if(iSumDelta <= NowPS.getSumDelta())
	    {
	      NowPS.setSumDelta( iSumDelta );
	      NowPS.setRowCoverage( makeRowGapCoverage(Row) );
	      NowPS.setColumnCoverage( ColumnCoverages.at(Col-1) );
	      NowPS.setBackDirection( 2 );
	      NowPS.setRowNumGap( prePS.getRowNumGap() );
	      NowPS.setColumnNumGap( prePS.getColumnNumGap()+1 );
	      
	      Nowcell.setTrueAndIndex(NumGap);
	    }
	}
      ++vit;
    }
}

void ShapeAlignment1::init()
{
  CELL<PathStatus> &firstcell=cells.at(0).at(0);
  firstcell.setTrue(0);
  firstcell.setIndex(0);
  PathStatus &PS=firstcell.getSS().at(0);
  PS.setSumDelta(0);

  int RowSize=cells.size();
  int ColSize=cells.at(0).size();
  for(int i=1; i<=min(MaxRowGap,RowSize-1); i++){ RowMove(i,0); }
  for(int j=1; j<=min(MaxColumnGap,ColSize-1); j++){ ColumnMove(0,j); }
  //for(int i=1; i<=MaxRowGap; i++){ RowMove(i,0); }
  //for(int j=1; j<=MaxColumnGap; j++){ ColumnMove(0,j); }
}

void ShapeAlignment1::Alignment()
{
  int RowSize=cells.size();
  int ColSize=cells.at(0).size();
  for(int i=1; i<RowSize; i++)
    {
      int fc=i-MaxRowGap;
      int ec=MaxColumnGap+i;
      for(int j=max(1,fc); j<=min(ColSize-1,ec); j++)
	{
	  RowMove(i,j);
	  ColumnMove(i,j);
	  DiagonalMove(i,j);
	}
    }
}

void ShapeAlignment1::TraceBack(vector<double> &Alignment_RowVector, vector<double> &Alignment_ColumnVector)
{
  CELL<PathStatus> *cell=&cells.back().back();
  sort( cell->getSS().begin(), cell->getSS().end() );

  PathStatus *PS=&cell->getSS().back();

  //cerr << PS->getSumDelta() << " : " << cell->getPSs().at(0).getSumDelta() << endl;

  Alignment_RowVector.push_back(PS->getRowCoverage());
  Alignment_ColumnVector.push_back(PS->getColumnCoverage());

  int Direction=PS->getBackDirection();
  int NumGap=PS->getNumGap();
  int i=cells.size()-1;
  int j=cells.at(0).size()-1;

  while(i != 0 or j != 0)
    {
      if(Direction == 0){i--; j--;}
      else if(Direction == 1){i--; NumGap-=1;}
      else if(Direction == 2){j--; NumGap-=1;}
      else{
	// cerr << "Direction value error : " << Direction << endl; 
	// cerr << i << " : " << j << " : " << NumGap << endl;
	// cerr << PS->getSumDelta() << endl;
	// exit(1);
        throw "Direction value error";
      }

      //cerr << i << " : " << j << " : " << NumGap << endl;

      cell=&cells.at(i).at(j);
      PS=&cell->getSS().at(NumGap);
      Alignment_RowVector.push_back( PS->getRowCoverage() );
      Alignment_ColumnVector.push_back( PS->getColumnCoverage() );
      Direction=PS->getBackDirection();
    }
}

float DPs::lnSum(const vector<float>& scores) const
{
  if(scores.size() == 1){return scores.at(0);}
  else
    {
      vector<float> scoresRemFront;
      scoresRemFront.assign(scores.begin()+1,scores.end());
      return scores.front()+log( 1+exp( lnSum(scoresRemFront)-scores.front() ) );
      }
}

float DPs::makefMscore(const vector<float>& scores, int Row, int Col) const
{
  return beta*GapScoreOfDiagonalMove(Row,Col) + lnSum( scores );
}

float DPs::makefRscore(const vector<float>& scores, int Row, int Col) const
{
  return beta*GapScoreOfRowMove(Row,Col) + lnSum( scores );
}

float DPs::makefCscore(const vector<float>& scores, int Row, int Col) const
{
  return beta*GapScoreOfColumnMove(Row,Col) + lnSum( scores );
} 

float DPs::makerscore(const vector<float> &scores) const
{
  return lnSum( scores );
} 

void KernelAlignment1::initForward()
{
  CELL<AlignmentKernelStatus> &firstcell=fcells.at(0).at(0);
  firstcell.setTrue(0);
  firstcell.setIndex(0);
  AlignmentKernelStatus &AKS=firstcell.getSS().at(0);
  AKS.setfMscore(0);
  //AKS.setfRscore(0);
  //AKS.setfCscore(0);
  AKS.setRowNumGap(0);
  AKS.setColumnNumGap(0);
  
  initForwardRow();
  initForwardCol();
}

void KernelAlignment1::initForwardRow()
{
  int RowSize=fcells.size();
  for(int i=1; i<=min(MaxRowGap,RowSize-1); i++)
    {
      CELL<AlignmentKernelStatus> &precell=fcells.at(i-1).at(0);
      CELL<AlignmentKernelStatus> &Nowcell=fcells.at(i).at(0);
      Nowcell.setTrue(i);
      Nowcell.setIndex(i);
      AlignmentKernelStatus &preAKS=precell.getSS().at(i-1);
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(i);
      NowAKS.setfRscore( makefRscore(preAKS.getfscores(),i,0) );
      //NowAKS.setfMscore(0);
      //NowAKS.setfCscore(0);
      NowAKS.setRowNumGap(i);
      NowAKS.setColumnNumGap(0);
    }
}

void KernelAlignment1::initForwardCol()
{
  int ColSize=fcells.at(0).size();
  for(int j=1; j<=min(MaxColumnGap,ColSize-1); j++)
    {
      CELL<AlignmentKernelStatus> &precell=fcells.at(0).at(j-1);
      CELL<AlignmentKernelStatus> &Nowcell=fcells.at(0).at(j);
      Nowcell.setTrue(j);
      Nowcell.setIndex(j);
      AlignmentKernelStatus &preAKS=precell.getSS().at(j-1);
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(j);
      NowAKS.setfCscore( makefCscore(preAKS.getfscores(),0,j) );
      //NowAKS.setfMscore(0);
      //NowAKS.setfRscore(0);
      NowAKS.setRowNumGap(0);
      NowAKS.setColumnNumGap(j);
    }
}

void KernelAlignment1::Forward()
{
  int RowSize=fcells.size();
  int ColSize=fcells.at(0).size();
  for(int i=1; i<RowSize; i++)
    {
      int fc=i-MaxRowGap;
      int ec=MaxColumnGap+i;
      for(int j=max(1,fc); j<=min(ColSize-1,ec); j++)
	{
	  CalfMscores(i,j);
	  CalfRscores(i,j);
	  CalfCscores(i,j);
	}
    }
  CalAllProb();
  if( isnan(AllProb) ){ cerr << "nan error in Forward" << endl; }
  if( isinf(AllProb) ){ cerr << "inf error in Forward" << endl; }
}

void KernelAlignment1::CalfMscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &precell=fcells.at(Row-1).at(Col-1);
  CELL<AlignmentKernelStatus> &Nowcell=fcells.at(Row).at(Col);

  vector<int>::const_iterator vit = precell.getIndex().begin();
  while( vit != precell.getIndex().end() ) 
    {
      AlignmentKernelStatus &preAKS=precell.getSS().at(*vit);
      int NumGap=preAKS.getNumGap();
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(NumGap);      

      Nowcell.setTrueAndIndex(NumGap);
      NowAKS.setfMscore( makefMscore(preAKS.getfscores(),Row,Col) );
      NowAKS.setRowNumGap( preAKS.getRowNumGap() );
      NowAKS.setColumnNumGap( preAKS.getColumnNumGap() );

      vit++;
    }
}

void KernelAlignment1::CalfRscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &precell=fcells.at(Row-1).at(Col);
  CELL<AlignmentKernelStatus> &Nowcell=fcells.at(Row).at(Col);

  vector<int>::const_iterator vit = precell.getIndex().begin();
  while( vit != precell.getIndex().end() ) 
    {
      AlignmentKernelStatus &preAKS=precell.getSS().at(*vit);
      int precellRowGap=preAKS.getRowNumGap();

      if( precellRowGap < MaxRowGap ) 
	{
	  int NumGap=preAKS.getNumGap()+1;
	  AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(NumGap);

	  Nowcell.setTrueAndIndex(NumGap);
	  NowAKS.setfRscore( makefRscore(preAKS.getfscores(),Row,Col) );
	  NowAKS.setRowNumGap( preAKS.getRowNumGap()+1 );
	  NowAKS.setColumnNumGap( preAKS.getColumnNumGap());
	}
      vit++;
    }
}

void KernelAlignment1::CalfCscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &precell=fcells.at(Row).at(Col-1);
  CELL<AlignmentKernelStatus> &Nowcell=fcells.at(Row).at(Col);

  vector<int>::const_iterator vit = precell.getIndex().begin();
  while( vit != precell.getIndex().end() ) 
    {
      AlignmentKernelStatus &preAKS=precell.getSS().at(*vit);
      int precellColumnGap=preAKS.getColumnNumGap();

      if( precellColumnGap < MaxColumnGap ) 
	{
	  int NumGap=preAKS.getNumGap()+1;
	  AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(NumGap);

	  Nowcell.setTrueAndIndex(NumGap);
	  NowAKS.setfCscore( makefCscore(preAKS.getfscores(),Row,Col) );
	  NowAKS.setRowNumGap( preAKS.getRowNumGap() );
	  NowAKS.setColumnNumGap( preAKS.getColumnNumGap()+1 );
	}
      vit++;
    }
}

void KernelAlignment1::initBackward()
{
  CELL<AlignmentKernelStatus> &firstcell=rcells.back().back();
  int MinGap = abs( (int)rcells.size() - (int)rcells.front().size() ); 

  for(int i=0; i<=MaxRowGap; i++)
    {
      firstcell.setTrueAndIndex(i*2+MinGap);
      AlignmentKernelStatus &AKS=firstcell.getSS().at(i*2+MinGap);
      AKS.setrMscore(0);
      AKS.setrRscore(0);
      AKS.setrCscore(0);
      AKS.setRowNumGap(i);
      AKS.setColumnNumGap(i+MinGap);
    }

  initBackwardRow();
  initBackwardCol();
}

void KernelAlignment1::initBackwardRow()
{
  int RowSize=rcells.size();
  int ColSize=rcells.back().size();
  for(int i=RowSize-2; i>=max(0,RowSize-MaxRowGap-1); i--)
    {
      CELL<AlignmentKernelStatus> &Nowcell=rcells.at(i).at(ColSize-1);
      CELL<AlignmentKernelStatus> &preRcell=rcells.at(i+1).at(ColSize-1);

      vector<int>::const_iterator vit = preRcell.getIndex().begin();
      while( vit != preRcell.getIndex().end() ) 
	{
	  AlignmentKernelStatus &preRAKS=preRcell.getSS().at(*vit);
	  int NumGap=preRAKS.getNumGap();
	  if(NumGap > 0)
	    {
	      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(NumGap-1);
	      Nowcell.setTrueAndIndex(NumGap-1);
	      
	      float prerRscore=preRAKS.getrRscore()+beta*GapScoreOfRowMove(i+1,ColSize-1);
	      NowAKS.setrMscore( prerRscore );
	      NowAKS.setrRscore( prerRscore );
	      NowAKS.setrCscore( prerRscore );
	      
	      NowAKS.setRowNumGap( preRAKS.getRowNumGap()-1 );
	      NowAKS.setColumnNumGap( preRAKS.getColumnNumGap() );
	      
	    }
	  vit++;
	}
    }
}

void KernelAlignment1::initBackwardCol()
{
  int RowSize=rcells.size();
  int ColSize=rcells.back().size();
  for(int j=ColSize-2; j>=max(0,ColSize-MaxColumnGap-1); j--)
    {
      CELL<AlignmentKernelStatus> &Nowcell=rcells.at(RowSize-1).at(j);
      CELL<AlignmentKernelStatus> &preCcell=rcells.at(RowSize-1).at(j+1);

      vector<int>::const_iterator vit = preCcell.getIndex().begin();
      while( vit != preCcell.getIndex().end() ) 
	{
	  AlignmentKernelStatus &preCAKS=preCcell.getSS().at(*vit);
	  int NumGap=preCAKS.getNumGap();
	  if(NumGap > 0)
	    {
	      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(NumGap-1);
	      Nowcell.setTrueAndIndex(NumGap-1);

	      float prerCscore=preCAKS.getrCscore()+beta*GapScoreOfColumnMove(RowSize-1,j+1);
	      NowAKS.setrMscore( prerCscore );
	      NowAKS.setrRscore( prerCscore );
	      NowAKS.setrCscore( prerCscore );

	      NowAKS.setRowNumGap( preCAKS.getRowNumGap() );
	      NowAKS.setColumnNumGap( preCAKS.getColumnNumGap()-1 );
	      
	    }
	  vit++;
	}
    }
}

void KernelAlignment1::Backward()
{
  int RowSize=rcells.size();
  int ColSize=rcells.back().size();
  for(int i=RowSize-2; i>=0; i--)
    {
      int fc=i-MaxRowGap;
      int ec=MaxColumnGap+i;
      for(int j=min(ColSize-2,ec); j>=max(0,fc); j--)
	{
	  Calrscores(i,j);
	}
    }

  float backprob=rcells.front().front().getSS().at(0).getrMscore();
  //cerr << "last backprob : " << backprob << endl;
  if( isnan(backprob) ){ cerr << "nan error in backward" << endl; }
  if( isinf(backprob) ){ cerr << "inf error in backward" << endl; }
}

void KernelAlignment1::Calrscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &Nowcell=rcells.at(Row).at(Col);
  CELL<AlignmentKernelStatus> &preMcell=rcells.at(Row+1).at(Col+1);

  vector<int>::const_iterator vit = preMcell.getIndex().begin();
  while( vit != preMcell.getIndex().end() ) 
    {
      AlignmentKernelStatus &preMAKS=preMcell.getSS().at(*vit);
      int NumGap=preMAKS.getNumGap();
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(NumGap);
      Nowcell.setTrueAndIndex(NumGap);

      vector<float> prerscores;
      setprerMscoreInVector(prerscores,Row,Col,*vit);
      setprerRscoreInVector(prerscores,Row,Col,*vit);
      setprerCscoreInVector(prerscores,Row,Col,*vit);
      if(prerscores.size() > 0)
	{
	  sort(prerscores.begin(),prerscores.end(),std::greater<float>());
	  float nowrscore=makerscore(prerscores);

	  NowAKS.setrMscore( nowrscore );
	  NowAKS.setrRscore( nowrscore );
	  NowAKS.setrCscore( nowrscore );
	}

      NowAKS.setRowNumGap( preMAKS.getRowNumGap() );
      NowAKS.setColumnNumGap( preMAKS.getColumnNumGap() );

      vit++;
    }
}

void KernelAlignment1::setprerMscoreInVector(vector<float> &prerscores, int Row, int Col, int index)
{
  CELL<AlignmentKernelStatus> &preMcell=rcells.at(Row+1).at(Col+1);
  AlignmentKernelStatus &preMAKS=preMcell.getSS().at(index);
  if(preMAKS.getrMb())
    {
      float prerMscore=preMAKS.getrMscore()+beta*GapScoreOfDiagonalMove(Row+1,Col+1);
      prerscores.push_back(prerMscore);
    }
}

void KernelAlignment1::setprerRscoreInVector(vector<float> &prerscores, int Row, int Col, int index)
{
  CELL<AlignmentKernelStatus> &preRcell=rcells.at(Row+1).at(Col);
  if( preRcell.getbool(index+1) )
    {
      AlignmentKernelStatus &preRAKS=preRcell.getSS().at(index+1);
      if(preRAKS.getrRb())
	{
	  float prerRscore=preRAKS.getrRscore()+beta*GapScoreOfRowMove(Row+1,Col);
	  prerscores.push_back(prerRscore);
	}
    }
}

void KernelAlignment1::setprerCscoreInVector(vector<float> &prerscores, int Row, int Col, int index)
{
  CELL<AlignmentKernelStatus> &preCcell=rcells.at(Row).at(Col+1);
  if( preCcell.getbool(index+1) )
    {
      AlignmentKernelStatus &preCAKS=preCcell.getSS().at(index+1);
      if(preCAKS.getrCb())
	{
	  float prerCscore=preCAKS.getrCscore()+beta*GapScoreOfColumnMove(Row,Col+1);
	  prerscores.push_back(prerCscore);
	}
    }

}

void KernelAlignment1::CalAllProb()
{
  vector<float> Lastfscores;
  CELL<AlignmentKernelStatus> &lastcell=fcells.back().back();
  
  for(int i=0; i<lastcell.getIndex().size(); i++)
    {
      int index=lastcell.getIndex().at(i);
      AlignmentKernelStatus &iAKS=lastcell.getSS().at(index);

      vector<float> fscores=iAKS.getfscores();
      for(int j=0; j<fscores.size(); j++){ Lastfscores.push_back(fscores.at(j)); }
    }
  sort(Lastfscores.begin(),Lastfscores.end(),std::greater<float>());

  AllProb = lnSum(Lastfscores);
  //cerr << "last forwardprog : " << AllProb << endl;
}

float KernelAlignment1::CalMpostProb(int Row, int Col)
{
  vector<float> Mscores;

  CELL<AlignmentKernelStatus> &fcell=fcells.at(Row).at(Col);
  CELL<AlignmentKernelStatus> &rcell=rcells.at(Row).at(Col);
  vector<int>::const_iterator vit = fcell.getIndex().begin();
  while( vit != fcell.getIndex().end() ) 
    {
      AlignmentKernelStatus &fAKS=fcell.getSS().at(*vit);
      AlignmentKernelStatus &rAKS=rcell.getSS().at(*vit);

      if(fAKS.getfMb() && rAKS.getrMb())
	{
	  Mscores.push_back( fAKS.getfMscore()+rAKS.getrMscore() );
	}
      vit++;
    }
  sort(Mscores.begin(),Mscores.end(),std::greater<float>());
  return lnSum(Mscores);
}

void KernelAlignment1::CalMPostProbs()
{
  int RowSize=fcells.size();
  int ColSize=fcells.at(0).size();
  
  for(int i=1; i<RowSize-1; i++)
    {
      int fc=i-MaxRowGap;
      int ec=MaxColumnGap+i;
      for(int j=max(1,fc); j<=min(ColSize-1,ec); j++)
	{
	  float MpostProb=CalMpostProb(i,j);
	  MpostProbs.at(i-1).at(j-1) = exp( MpostProb-AllProb );
	  //cerr << "MpostProb : " << MpostProb << " : AllProb : " << AllProb << endl;
	  //cerr << i << " : " << j << " Prob : " <<  MpostProbs.at(i-1).at(j-1) << endl;
	  //if( MpostProbs.at(i-1).at(j-1) > 1.0){cerr << "error in Calculate AlignmentProbs. if beta >= 0, beta should be < 0." << endl; exit(1);}
	  if( MpostProbs.at(i-1).at(j-1) > 1.0){ MpostProbs.at(i-1).at(j-1) =1.0; }
	}
    }
}

void KernelAlignment2::initForward() 
{
  CELL<AlignmentKernelStatus> &firstcell=fcells.at(0).at(0);
  AlignmentKernelStatus &AKS=firstcell.getSS().at(0);
  AKS.setfMscore(0);
  //AKS.setfRscore(0);
  //AKS.setfCscore(0);
  
  initForwardRow();
  initForwardCol();
}

void KernelAlignment2::initForwardRow()
{
  int RowSize=fcells.size();
  for(int i=1; i<=RowSize-1; i++)
    {
      CELL<AlignmentKernelStatus> &precell=fcells.at(i-1).at(0);
      CELL<AlignmentKernelStatus> &Nowcell=fcells.at(i).at(0);
      AlignmentKernelStatus &preAKS=precell.getSS().at(0);
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);
      NowAKS.setfRscore( makefRscore(preAKS.getfscores(),i,0) );
      //NowAKS.setfMscore(0);
      //NowAKS.setfCscore(0);
    }
} 

void KernelAlignment2::initForwardCol()
{
  int ColSize=fcells.at(0).size();
  for(int j=1; j<=ColSize-1; j++)
    {
      CELL<AlignmentKernelStatus> &precell=fcells.at(0).at(j-1);
      CELL<AlignmentKernelStatus> &Nowcell=fcells.at(0).at(j);
      AlignmentKernelStatus &preAKS=precell.getSS().at(0);
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);
      NowAKS.setfCscore( makefCscore(preAKS.getfscores(),0,j) );
      //NowAKS.setfMscore(0);
      //NowAKS.setfRscore(0);
    }
}

void KernelAlignment2::Forward()
{
  int RowSize=fcells.size();
  int ColSize=fcells.at(0).size();
  for(int i=1; i<RowSize; i++)
    {
      for(int j=1; j<=ColSize-1; j++)
	{
	  CalfMscores(i,j);
	  CalfRscores(i,j);
	  CalfCscores(i,j);
	}
    }
  CalAllProb();
  if( isnan(AllProb) ){ cerr << "nan error in Forward" << endl; }
  if( isinf(AllProb) ){ cerr << "inf error in Forward" << endl; }
}

void KernelAlignment2::CalfMscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &precell=fcells.at(Row-1).at(Col-1);
  CELL<AlignmentKernelStatus> &Nowcell=fcells.at(Row).at(Col);

  AlignmentKernelStatus &preAKS=precell.getSS().at(0);
  AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);      
  
  NowAKS.setfMscore( makefMscore(preAKS.getfscores(),Row,Col) );
}

void KernelAlignment2::CalfRscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &precell=fcells.at(Row-1).at(Col);
  CELL<AlignmentKernelStatus> &Nowcell=fcells.at(Row).at(Col);

  AlignmentKernelStatus &preAKS=precell.getSS().at(0);
  AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);
  
  NowAKS.setfRscore( makefRscore(preAKS.getfscores(),Row,Col) );
}

void KernelAlignment2::CalfCscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &precell=fcells.at(Row).at(Col-1);
  CELL<AlignmentKernelStatus> &Nowcell=fcells.at(Row).at(Col);

  AlignmentKernelStatus &preAKS=precell.getSS().at(0);
  AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);

  NowAKS.setfCscore( makefCscore(preAKS.getfscores(),Row,Col) );
}

void KernelAlignment2::initBackward()
{
  CELL<AlignmentKernelStatus> &firstcell=rcells.back().back();
  AlignmentKernelStatus &AKS=firstcell.getSS().at(0);
  AKS.setrMscore(0);
  AKS.setrRscore(0);
  AKS.setrCscore(0);

  initBackwardRow();
  initBackwardCol();
}

void KernelAlignment2::initBackwardRow()
{
  int RowSize=rcells.size();
  int ColSize=rcells.back().size();
  for(int i=RowSize-2; i>=0; i--)
    {
      CELL<AlignmentKernelStatus> &Nowcell=rcells.at(i).at(ColSize-1);
      CELL<AlignmentKernelStatus> &preRcell=rcells.at(i+1).at(ColSize-1);
      
      AlignmentKernelStatus &preRAKS=preRcell.getSS().at(0);
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);
	      
      float prerRscore=preRAKS.getrRscore()+beta*GapScoreOfRowMove(i+1,ColSize-1);
      NowAKS.setrMscore( prerRscore );
      NowAKS.setrRscore( prerRscore );
      NowAKS.setrCscore( prerRscore );      
    }
}

void KernelAlignment2::initBackwardCol()
{
  int RowSize=rcells.size();
  int ColSize=rcells.back().size();
  for(int j=ColSize-2; j>=0; j--)
    {
      CELL<AlignmentKernelStatus> &Nowcell=rcells.at(RowSize-1).at(j);
      CELL<AlignmentKernelStatus> &preCcell=rcells.at(RowSize-1).at(j+1);

      AlignmentKernelStatus &preCAKS=preCcell.getSS().at(0);
      AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);

      float prerCscore=preCAKS.getrCscore()+beta*GapScoreOfColumnMove(RowSize-1,j+1);
      NowAKS.setrMscore( prerCscore );
      NowAKS.setrRscore( prerCscore );
      NowAKS.setrCscore( prerCscore );      
    }
}

void KernelAlignment2::Backward()
{
  int RowSize=rcells.size();
  int ColSize=rcells.back().size();
  for(int i=RowSize-2; i>=0; i--)
    {
      for(int j=ColSize-2; j>=0; j--)
	{
	  Calrscores(i,j);
	}
    }

  float backprob=rcells.front().front().getSS().at(0).getrMscore();
  //cerr << "last backprob : " << backprob << endl;
  if( isnan(backprob) ){ cerr << "nan error in backward" << endl; }
  if( isinf(backprob) ){ cerr << "inf error in backward" << endl; }
}

void KernelAlignment2::Calrscores(int Row, int Col)
{
  CELL<AlignmentKernelStatus> &Nowcell=rcells.at(Row).at(Col);
  CELL<AlignmentKernelStatus> &preMcell=rcells.at(Row+1).at(Col+1);

  AlignmentKernelStatus &preMAKS=preMcell.getSS().at(0);
  AlignmentKernelStatus &NowAKS=Nowcell.getSS().at(0);
  
  vector<float> prerscores;
  setprerMscoreInVector(prerscores,Row,Col);
  setprerRscoreInVector(prerscores,Row,Col);
  setprerCscoreInVector(prerscores,Row,Col);

  sort(prerscores.begin(),prerscores.end(),std::greater<float>());
  float nowrscore=makerscore(prerscores);

  NowAKS.setrMscore( nowrscore );
  NowAKS.setrRscore( nowrscore );
  NowAKS.setrCscore( nowrscore );
}

void KernelAlignment2::setprerMscoreInVector(vector<float> &prerscores, int Row, int Col)
{
  CELL<AlignmentKernelStatus> &preMcell=rcells.at(Row+1).at(Col+1);
  AlignmentKernelStatus &preMAKS=preMcell.getSS().at(0);
  float prerMscore=preMAKS.getrMscore()+beta*GapScoreOfDiagonalMove(Row+1,Col+1);
  prerscores.push_back(prerMscore);
}

void KernelAlignment2::setprerRscoreInVector(vector<float> &prerscores, int Row, int Col)
{
  CELL<AlignmentKernelStatus> &preRcell=rcells.at(Row+1).at(Col);
  AlignmentKernelStatus &preRAKS=preRcell.getSS().at(0);
  float prerRscore=preRAKS.getrRscore()+beta*GapScoreOfRowMove(Row+1,Col);
  prerscores.push_back(prerRscore);
}

void KernelAlignment2::setprerCscoreInVector(vector<float> &prerscores, int Row, int Col)
{
  CELL<AlignmentKernelStatus> &preCcell=rcells.at(Row).at(Col+1);
  AlignmentKernelStatus &preCAKS=preCcell.getSS().at(0);
  float prerCscore=preCAKS.getrCscore()+beta*GapScoreOfColumnMove(Row,Col+1);
  prerscores.push_back(prerCscore);
}

void KernelAlignment2::CalAllProb()
{
  CELL<AlignmentKernelStatus> &lastcell=fcells.back().back();  
  AlignmentKernelStatus &iAKS=lastcell.getSS().at(0);

  vector<float> Lastfscores=iAKS.getfscores();
  sort(Lastfscores.begin(),Lastfscores.end(),std::greater<float>());

  AllProb = lnSum(Lastfscores);
  //cerr << "last forwardprog : " << AllProb << endl;
}

float KernelAlignment2::CalMpostProb(int Row, int Col)
{
  vector<float> Mscores;

  CELL<AlignmentKernelStatus> &fcell=fcells.at(Row).at(Col);
  CELL<AlignmentKernelStatus> &rcell=rcells.at(Row).at(Col);
  AlignmentKernelStatus &fAKS=fcell.getSS().at(0);
  AlignmentKernelStatus &rAKS=rcell.getSS().at(0);

  if(fAKS.getfMb() && rAKS.getrMb())
    {
      Mscores.push_back( fAKS.getfMscore()+rAKS.getrMscore() );
    }
  sort(Mscores.begin(),Mscores.end(),std::greater<float>());
  return lnSum(Mscores);
}

void KernelAlignment2::CalMPostProbs()
{
  int RowSize=fcells.size();
  int ColSize=fcells.at(0).size();
  
  for(int i=1; i<RowSize-1; i++)
    {
      for(int j=1; j<=ColSize-1; j++)
	{
	  float MpostProb=CalMpostProb(i,j);
	  MpostProbs.at(i-1).at(j-1) = exp( MpostProb-AllProb );
	  //cerr << "MpostProb : " << MpostProb << " : AllProb : " << AllProb << endl;
	  //cerr << i << " : " << j << " Prob : " <<  MpostProbs.at(i-1).at(j-1) << endl;
	  //if( MpostProbs.at(i-1).at(j-1) > 1.0){cerr << "error in Calculate AlignmentProbs. if beta >= 0, beta should be < 0." << endl; exit(1);}
	  if( MpostProbs.at(i-1).at(j-1) > 1.0){ MpostProbs.at(i-1).at(j-1) =1.0; }
	}
    }
}

////////////////////////////////////////////////////////////////////////////////////////

void
normalize_vector(const std::vector<double>& v, std::vector<double>& w)
{
  double max_v = *std::max_element(v.begin(), v.end());
  w=v;
  for (std::vector<double>::iterator x=w.begin(); x!=w.end(); ++x)
    *x /= max_v;
}

void
ShapeProbAlign::
calculate_posterior(const std::vector<double>& cov1, const std::vector<double>& cov2,
                    std::vector<std::vector<float> >& posterior) const
{
  //  if (cov1.size()>=cov2.size())
  if (cov1.size()>cov2.size())
  {
    std::vector<std::vector<float> > posterior_temp;
    calculate_posterior(cov2, cov1, posterior_temp);
    posterior.resize(posterior_temp[0].size(), std::vector<float>(posterior_temp.size(), 0.0));
    for (unsigned int i=0; i!=posterior_temp.size(); ++i)
      for (unsigned int j=0; j!=posterior_temp[0].size(); ++j) //i!= - j!=
        posterior[j][i]=posterior_temp[i][j];
  }
  else
  {
    std::vector<double> v1, v2;
    normalize_vector(cov1, v1);
    normalize_vector(cov2, v2);
    int max_gap;
    if (max_gap_<0)             // "flow" mode
      max_gap = (int) v2.size()*alpha_;
    else                        // "fix" model
      max_gap = max_gap_;

    if(alpha_ < 1.0)
      {
	KernelAlignment1 AK1(v1, v2, max_gap, beta_, both_end_zero_); //reverse beta_ both_end_zero_
	AK1.initForward();
	AK1.Forward();
	AK1.initBackward();
	AK1.Backward();
	AK1.CalMPostProbs();
	posterior=AK1.getMpostProbs();
      }
    else // alpha_ == 1.0
      {
	KernelAlignment2 AK2(v1, v2, 1, beta_, both_end_zero_); //reverse beta_ both_end_zero_
	AK2.initForward();
	AK2.Forward();
	AK2.initBackward();
	AK2.Backward();
	AK2.CalMPostProbs();
	posterior=AK2.getMpostProbs();
      }
  }
}

void
ShapeProbAlign::
calculate(const std::string& seq1, const std::string& seq2, MP& mp)
{
  throw "not supported";
}

void
ShapeProbAlign::
calculate(const std::vector<Fasta>& fa, std::vector<std::vector<MP> >& mp)
{
  const uint N=cov_.size();
  if (cov_.size()!=fa.size())
  {
    throw "Error: mismatch the number of coverage vectors";
    return;
  }
  for (uint i=0; i!=N; ++i)
  {
    if (cov_[i].size()!=fa[i].size())
    {
      std::ostringstream s;
      s << "Error: mismatch of the length of " << i << "-th sequence and cov vector: "
        << fa[i].size() << "!=" << cov_[i].size();
      throw s.str().c_str();
      return;
    }
  }

  mp.resize(N, std::vector<MP>(N));
  for (uint i=0; i!=N; ++i)
  {
    mp[i][i].resize(cov_[i].size());
    for (uint x=0; x!=cov_[i].size(); ++x)
      mp[i][i][x].push_back(std::make_pair(x, 1.0f));
    for (uint j=i+1; j!=N; ++j)
    {
      mp[i][j].resize(cov_[i].size());
      std::vector<std::vector<float> > posterior;
      calculate_posterior(cov_[i], cov_[j], posterior);
      for (uint x=0; x!=cov_[i].size(); ++x)
      {
        for (uint y=0; y!=cov_[j].size(); ++y)
        {
          if (posterior[x][y]>threshold())
            mp[i][j][x].push_back(std::make_pair(y, posterior[x][y]));
        }
      }
    }
  }
}
